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Radial stability of a family of anisotropic Hernquist models 
with and without a supermassive black hole 
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ABSTRACT 

We present a method to investigate the radial stability of a spherical anisotropic 
system that hosts a central supermassive black hole (SBH). Such systems have never 
been tested before for stability, although high anisotropies have been considered in 
the dynamical models that were used to estimate the masses of the central putative 
supermassive black holes. A family of analytical anisotropic spherical Hernquist models 
with and without a black hole were investigated by means of .ZV-body simulations. A 
clear trend emerges that the supermassive black hole has a significant effect on the 
overall stability of the system, i.e. an SBH with a mass of a few percent of the total 
mass of the galaxy can prevent or reduce the bar instabilities in anisotropic systems. 
Its mass not only determines the strength of the instability reduction, but also the time 
in which this occurs. These effects are most significant for models with strong radial 
anisotropies. Furthermore, our analysis shows that unstable systems with similar SBH 
but with different anisotropy radii evolve differently: highly radial systems become 
oblate, while more isotropic models tend to form into prolate structures. In addition 
to this study, we also present a Monte-Carlo algorithm to generate particles in spherical 
anisotropic systems. 

Key words: Stellar dynamics - methods : N-body simulations - galaxies : kinematics 
and dynamics 



1 INTRODUCTION 



Nowadays it is accepted that almost every galaxy hosts a 
central supermassive black hole (SBH) at its core. Since 
the kinematical discovery of the first SBH with the Hub- 
ble Space Telescope (HST), extensive studies have been 
carried out by many groups that investigate the demog- 
raphy of SBHs and the effect of the SBHs on their en- 
vironment. The most popular discoveries are the correla- 
tions between the mass of the SBH (Mbh) and respec- 
tively the total blue magnit ude Lb of the hot stellar com - 
ponent in which it resides (|Kormendv fc Richstone 19951 ). 
the central velocity disper s ion of the hot stellar comp onent 
jFerrarese fc Merritt 200(il ; iGebhardt et al. 2000h . the cen- 
tral lig ht concentration C( a) or equivalent the Sersic in- 
dex n l|Graham et al. 2001) and the maximum rotational 
velocity of the galaxy jFerrarese 20021 ; iBaes et al. 20031 ; 
iPizzella et al. 20051 ; iBuvle et al. 200rj ). These relations have 
been calibrated with the known masses of the SBHs of the 
nearest galaxies, that mostly have been derived by means of 
either stellar or gas kinematics. 

Sophisticated axisymmetric 3-integral dynamical mod- 
els that allow a variation in mass-to-light ratio and 
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anisotropy as a function of radius have been obtained by fits 
to the line-of-sight velocity distributions (LOSVDs) in the 
galaxies, which were derived primarily from high-resolution 
spectra taken with the HST 



LOSVD (x,y,v z 



MI 



F(r, v) dz dv x dv y , 



(1) 



where F(r,v) denotes the stellar distribution function (DF) 
and p p stands for the projected mass density at position 
(x,y). The accuracy of the applied dynamical models to the 
observed stellar kinematics is still improving steadily and is 
reflected on the complexity of the DFs. Despite this posi- 
tive progress on the dynamical front, very few anisotropic 
dynamical models o f a galactic nucleus hav e been tested for 
dynamical stability jFerrarese fc Ford 20051 ). One of the rea- 
sons for this is the complexity of the distribution functions, 
which are mostly numerically derived. Hence, to simulate 
these numerical DFs one normally approximates numerically 
the solution of the Jeans equations to derive the velocity 
dispersion profile and then uses Gaussians to provide local 
velocity distributions. It is now known from recent simu- 
lations of galactic s ystems that this metho d causes serious 
numerical artifacts (Ka zantzidis et al. 2004h . 

So far the only theoretical analytical systems that con- 
tain a SBH are derived by Ciotti (1996) and Baes et al. 
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(2004, 2005) where the attention is drawn primarily to the 
Hernquist model since this is the best-known approximation 
to the Sersic profiles that are observed in bulges and ellipti- 
cal galaxies, and by Stiavelli (1998) where the distribution 
function of a stellar system around an SBH is derived from 
statistical mechanic considerations. Ciotti (1996) initially 
starts with a 2-component system containing the luminous 
and dark matter and creates both isotropic and anisotropic 
(based on the Osipkov-Merritt strategy) systems. The dark 
matter halo (also represented by a Hernquist model) can be 
transformed into a central SBH by setting the core radius 
to zero. 

In this article we present for the first time the re- 
sults of a dynamical stability investigation of spherical sys- 
tems containing an SBH, as a function of the mass of the 
SBH and the anisotropy radius of the system. We fill fo- 
cus primarily on the so-called radial orbit instability of 
radially anisotropic, spherically s ymmetric stellar systems 
l|Henon 19731 ; ICincotta et al. 19 96). A beautiful mathemati- 
cal and physical explanation for this instability can be found 
in Palmer (1994) and Merritt (1987) and references therein. 
As a bar grows, it saps angular momentum from the stars 
as they precess through the bar. Stars on low-angular mo- 
mentum orbits are trapped into resonance and strengthen 
the bar, making it possible to also trap stars with higher 
angular momentum into resonance and so on. Eventually, 
the triaxial force field of the bar becomes dominant and the 
initially rosette-shaped orbits in a spherical potential are 
transformed into box orbits. Radial orbits and box orbits 
bring stars close to the center of the galaxy where they can 
be diverted from their orbits by the spherically symmetric 
force field of a central massive black hole. This in turn may 
weaken the bar over time, or even prevent it from forming 
in the first place, clearly proving the relevance of a study 
such as the one presented here. 

In Section 2 we describe a Monte-Carlo algorithm that 
we developed to generate the initial conditions for the mod- 
els, together with our iV-body code and technique for inves- 
tigating the stability. We present in Section 3 the results of 
a stability investigation of a family of anisotropic Hernquist 
mo dels without an SBH, w ith different anisotropy behav- 
ior (|Baes fc Deionghe 20021 ). In Section 4 we describe the 
Osipkov-Merritt Hernquist models with a central SBH, in- 
troduced by Ciotti (1996). We investigate the stability of 
these systems in detail in Section 5, comparing them with 
the according models without an SBH. We perform this in a 
2-parameter space as a function of the anisotropy radius r a 
and the mass of the central SBH fi. In Section 6 we present 
our final results and conclusions. 



2 COMPUTATIONAL METHOD 

2.1 Definition of the Hernquist models 

First of all we introduce some general characteristics of the 
models in our dynamical study. All systems are based on th e 
spherical Hernquist potential-density pair dHernquist 1990T ). 
including a supermassive black hole in the centre. Given this 
mass profile, we shall investigate several distribution func- 
tions (DFs) consistent with the density outside the center, 
and which we will refer to as the stellar component. If we 




Figure 1. Visualization of an isotropic Hernquist system with 
fi = 0.1 and our approximation with cells. After 8 subdivisions 
991 cells were constructed, with a total phase-space volume of 
Vg = 1.533, while the real total stellar DF volume is 1. Thus, 
the ratio of rejected to accepted particles is 0.533 and on average 
~ 35% of all randomly chosen test particles in the cell volume will 
be rejected, resulting in a highly efficient Monte-Carlo simulation. 



denote the total stellar mass by M s , we can write the total 
mass as M to t ~ M s (l + fi), where the fractional quantity 
determines the SBH mass /iM s . In our subsequent analysis 
we will work in dimensionless units G = M s = 1, so that the 
gravitating binding potential and the density are given by 



i/)(r) = 
p(r) = 



1 



r r 



27rr(l + r) 3 



(r > 0). 



(2) 
(3) 



We will also express the time-steps in our A-body code 
(the time between two successive calculations) in dimension- 
less units of half-m ass dynamical time, which is defined as 
the dynamical time l|Binnev fc Tremaine 19871 ) at the stellar 
half-mass radius: 



3tt 
16Gp' 



where 

__ 3Af(n /2 ) 

P 4 ^l/2 



(4) 



(5) 



For a Hernquist model with fi — the half-mass dynamical 
time and the half-mass radius are 



n/2 



V2 

7 

2 

1 + V2 



1 + V2 



3/2 



(6) 
(7) 



We will also use these units for models with an SBH. 

A conversion to observational units can be obtained 
through the close sim ilarity between the Hernquist and De 
Vaucouleurs profiles (|Hernquist 199Gn . with r 1 / 2 ~ 1.33r e 
where r e is the effective radius. Then a physical length, time 
and velocity are found by the scaling relations 



r u r, 



(8) 
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Figure 2. The figures show the relevant parameters (p(r), o>(r), cxt(r) and /3(r)) of the outcome of a Monte Carlo sampling of a 
Hernquist system with fi = 0.1 and r a = 1. The continuous lines denote the theoretical model, the discrete data represent 10 5 simulated 
particles, binned and with error-bars. 



n 



GM, 



GMu 



■Th, 



-v, 



(9) 



(10) 



with M a the stellar mass, M to t the total mass and f u = 
(1.33r e )/(l + %/2), expressed in physical units. 

Every model was simulated by means of 10 5 equal -mass 
particles that all follow the distribution function of the sys- 
tem and are contained within a sphere of radius rt = 2000 
which encloses about 99.9% of the stellar mass of the sys- 
tem. We performed the simulations for 50 dynamical times, 
and used the values of the axis ratios c/a and b/a during 
this time as (in)stability indicators (see Section [2. 41) . 



2.2 Constructing the data sets 

Since we will investigate our models by means of iV-body 
simulations, the first objective is to obtain representative 
discrete data sets from the considered distribution func- 
tions. Each of these functions describes a spherical mass 
distribution in a dynamical system governed by a gravita- 
tional binding potential ip{r) > 0, which implies that they 
can be expressed as functions F(E, L) of the binding en- 



ergy E — V>( r ) — \ v "r ~ \ v t an d the modulus of the angular 
momentum L = rv t , with v r and Vt the radial and tangen- 
tial velocity components, respectively. Isotropic DFs can be 
reduced to F(E). 

In order to extract discrete data samples from the dis- 
tributions, we need to simulate random particles uniformly 
in the phase-space enclosed by the DFs. To this aim we used 
a Monte-Carlo simulator, developed by one of the authors 
(E.V.H.). The procedure works as follows: we write each DF 
as F(r,v r ,vt) and we consider a 4-dimensional grid space 
with (r, v r ,Vt) as abscissae and the function values on the 
ordinate axis. 

We start with a single cell in this space, extending from 
the origin to a boundary (rj,, Wr,6i v t,t) (where rt is chosen to 
be sufficiently large, and v r ,b = ft, 6 = \/2ip(0) ), and with 
the ordinate set at the (known or estimated) DF maximum 
fb- These boundaries (for infinite values a sufficiently large 
value is taken, see further) enclose a 7-dimensional phase- 
space volume 



Vi 



4tt 3 



2v r ,i 



(11) 



In the second step we attempt to split the cell into 
8 sub-cells with different ordinates (i.e. the up to that 
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point known function maxima in each cell) . Therefore a co- 
ordinate (r s , v rtS , Vt, s ) is sought to serve as the common cor- 
ner point in the abscissae for these sub-cells: starting in the 
cell center, the total phase-space volume of the originating 
sub-cells is calculated, and through a number of iterations 
the cell is scanned for a better splitting point, i.e. which 
minimizes this volume. In this manner, the original cell is 
being split as efficiently as possible into 8 new cells, adding 
up to a new total volume 

8 
i=l 

8 8tt 2 

= '^^{rl- l -rl. i ){v rib . l -v r ^ a , i )(v 2 t , b . i -vl ia . i )f b . i , (12) 
i=i 6 

which is a better approximation to the real DF volume. Here, 
for a cell i we denoted Vi-$ its volume, (r a -i,v r ,a;i,Vt,a;i) and 
(fb;i,Vr,b;i, v t,b;i) its lower and upper bounds in the abscissae, 
and fb-i its maximum DF value. 

Next, each cell in our grid is examined according to the 
procedure above and split if it leads to a significant decrease 
in the total volume. Thus, after the examination of every 
cell, a new volume V3 is obtained. This loop is repeated un- 
til after M steps the phase-space volume Vm has converged 
sufficiently close to the real volume. Typically, in our sim- 
ulations, the cells cover a volume that is a factor 1.5 to 5 
larger than the model's actual phase-space volume; a fur- 
ther refinement is unnecessary, since constructing more cells 
would be more time-consuming than actually generating our 
desired number (10 5 ) of particles (see below). If the grid is 
successfully constructed, F(r,v r ,Vt) is entirely enveloped by 
a set of 4-dimensional grid cells. 

Now we can proceed to a classical rejection Monte-Carlo 
(MC) simulation (in the remainder, we refer to setting up 
the initial conditions of a DF as an "MC simulation"). To 
generate a data point n, first a value V n is randomly chosen 
between and Vm . We can associate this value with a unique 
cell j and an ordinate f n for which 

J2 V M;i V M ;i, (13) 

i=l i=l 

and 

j-i 

Vn = Vm <* + 

i=l 

-^(rij-rl^)(v r fi.j-v r>a .j)(vt >b . i j-vt, a -,j)f n - (14) 

Then, in cell j the co-ordinates r 3 a .j < r 3 < rf ;j , w r ,a;j ^ 
v r ,n < v r ,b;j and v\ a .j < Vt : „ < v t,b;j are randomly gen- 
erated. Thus, a point (r n , tv,n> Vt,n, fn) is uniformly cho- 
sen in the 7-dimensional phase-space volume Vm- Now, if 
fn < F(r n ,v r ,n,v t ,n), the co-ordinate (r», V r , n , «t,n) is ac- 
cepted as a valid data point, otherwise it is rejected. Fur- 
thermore, if fb-j < F(r n ,v rtn ,vt,n), the cell volume is ac- 
cordingly increased to the new maximum, so the grid keeps 
being improved. 

In this manner we construct a data set of TV accepted 
co-ordinates inside the chosen radius r b which follow the 
distribution. The MC simulation is regarded successful if the 
cell volumes have changed negligibly (if the relative change 
of the total volume is smaller than 10 -3 ) during the MC 



simulation. If not, a new MC simulation with the final grid 
(with volume Vm+n) is necessary. Also, if the ratio between 
rejected and accepted points is very large, causing the MC 
simulation to be slow, the grid might have to be refined 
further (as aforementioned, we stop refining the grid once 
the cells cover a volume that is a factor 1.5 to 5 larger than 
the model's actual phase-space volume). 

Finally, every co-ordinate (r n , v r , n , v t,n) has to be con- 
verted into a phase-space point (x n ,yn, z n , v x ,„, %,n, 
This is done by uniformly simulating the surface of a sphere 
with radius r n (creating (x n ,y„, z n )), a circle with radius 
vt,n (creating («e, n ,w^,n)) and the sign of v r , n - The veloci- 
ties can then be transformed into the appropriate Cartesian 
co-ordinates. For isotropic functions F(E) the grid abscis- 
sae simplify to the 2-dimensional (r, v) space, and the entire 
procedure is analogous. 

Our method has several advantages: the construction 
of a grid and the subsequent MC simulation of points 
is straightforward, fast, accurate and generally applicable. 
This contrasts with algorithms that require integrations and 
inversions of DFs, which can experience numerical problems 
with intricate functions. Also, no intermediate steps are re- 
quired (e.g. simulating the density first before assigning ve- 
locities to each particle) and once a grid is made for a model, 
it can be re-used to generate an arbitrary number of parti- 
cles. Moreover, since a peak can be adequately isolated by 
a cell, infinite ranges in the co-ordinate space or the DF 
values can be approximated by choosing appropriate large 
boundary values. 

As an example we show in Fig. [T] the constructed cells 
for an isotropic Hernquist system with a central SBH of p = 
0.1 (see Section 4). A simulated data sample (10 5 accepted 
particles) for an anisotropic Hernquist system with a central 
SBH of p — 0.1 and an anisotropy radius of r a — 1 is shown 
in Fig. [2] In all our MC simulations, we truncate the infinite 
boundary radius at r b = 2000. For the DFs with an infinite 
maximum, we set f b = 10 15 , and for the SBH-models we 
set the maximum velocity at the arbitrarily large value v b = 
10 15 (these values are in fact much larger than needed. In 
reality, no particle is ever assigned such a high DF value or 
initial velocity and never reaches such high velocities during 
the subsequent TV-body simulations). 



2.3 TV-body code 

We studied the stability of our models by using an TV-body 
code that is based on the "self-consistent field" method 
l|Hernquist fc Ostriker 19921 ). This method relies on the se- 
ries expansion in a bi-orthogonal spherical basis set for both 
the density and gravitational potential 

p(r,6,4>) = ^2A nim p n lm{r,9,4>) 

nlm 

= J2 Anlm Pm{r)Yim{9,4>), (15) 

nlm 

$>(r,6,4>) = ^4i ra $ nim (r,6»,0) 

nlm 

= ^Ani m ^ni{r)Y lm {e,cj)), (16) 

nlm 
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Figure 3. (a) The initial particle positions of a Hernquist model with an increasing anisotropy with /3 = 0.5 and A = 5. 80% of the 
total mass is shown in the figure, (b) The density distribution after 10 half-mass dynamical times. A bar is clearly visible, (c) The axis 
ratio c/a plotted as a function of time. As can be seen from both the density distributions and the axis ratio, an elliptical bar is created 
indicating that the system is unstable. 



where Yi m (6, cj>) are the spherical harmonics. Some freedom 
is considered for this expansion since (p n ;(r), $„;(?••)) can 
have different forms (e.g. Plummer model, Bessel functions, 
spherical harmonic functions), however here we will use a 
form similar to the Hernquist model due to its trivial con- 
nection with our anisotropic systems that we wish to exam- 
ine: 



Pm(r) 



(2!+3/2) 



(0, 



V^r(l+ r ) 2 '+ 3 " 
-2v ^ , r _l a+1 Cf +3/2) (fl, 



(1 + r) 



(17) 
(18) 



where K n i is a normalization constant, £ = (r— l)/(r+l) and 
q(21+3/2) ^ are Q e g enDauer polynomials (e.g. Szego (1939), 
Sommerfeld (1964)). The coefficients A n i m can be calculated 
by means of all the particles that describe the DF of our 
system (see Hernquist & Ostriker (1992) for more details). 
The spherical accelerations for each particle are found by 
taking the gradient of the potential (eq. 1 16p . Finally new 
positions and velocities are derived with the use of an in- 
tegrator which is equivalent to t he standard tim e-centred 
leapfrog j Allen fc Tildeslev 19921 ; but et al. lmEt) , 

1 2 

Xi+i = Xi + Atvi + -At en, (19) 
Vi+i = Vi + ^At(a; + Oi+x). (20) 

The indices n, (I, m = —I. ..I) are indirectly an indication 
for the accuracy of the simulation for respectively radial 
and tangential motion, since they determine the number of 
terms in the expansion (see Section 5.2 Hernquist & Ostriker 
(1992) for a statistical analysis). For the systems without the 
SBH we find that n max =4 and lmax=i assures a total en- 
ergy conservation of better than ~ 1CP 6 over 50 half-mass 
dynamical times Th and still allows a low CPU time per 
A^-body time-step (the time between two successive calcu- 
lations) of At — Th/416 w 0.02. For the systems with an 
SBH, we used n ma x = 6, l ma x = 2 when fi < 0.05, and 
rimax = 8, l m ax = 4 for larger values of fi. The gravitational 
effect of the SBH is added analytically by an extra radial 
acceleration proportional to the mass of the SBH. To avoid 
numerical divergences when particles pass close to the SBH, 



we included a softening to this acceleration, i.e. — n/(r 2 + e 2 ) 
with softening length e = 0.05. At this radius the dynamical 
crossing time of a particle is Th = 0.37, which is still no- 
tably larger than our time-step of 0.02. Adding this softening 
causes a discrepancy in the treatment of the SBH potential 
between the analytical models and the N-body code. As a 
consequence the initial conditions are not exactly in dynam- 
ical equilibrium so that the systems develop transient radial 
motions to adjust their density profile. However, this effect 
is marginal and does not influence the overall results. In all 
simulations with an SBH the energy is conserved better than 
1% over 50 half-mass dynamical times. 

In order to check the robustness of our results, we per- 
formed two kinds of tests. We (i) re-ran a number of sim- 
ulations with different, smaller time-steps, and (%%) we per- 
formed simulations with higher n max and lmax values. A 
detailed comparison of these extra runs with the original 
simulations shows that our results and conclusions do not 
change : the variation of the global instability indicators, 
such as axis ratios or 2K r /Kt, as a function of time are 
essentially the same. 



2.4 Quantifying the instabilities 

When a system is unstable, it tends to create a bar fea- 
ture at its center (see Fig. [3J) which rough ly has an el- 
lipsoidal shape. As noted by other authors |Merritt 19871 ; 
iPalmer fc Papaloizou 1987h . the physical cause of instabil- 
ity is similar to th at of the formation of a bar in a disc 
l|Lvnden-Bell 1973 ) , where a small perturbation changes the 
orbits with a lower angular momentum (initially process- 
ing ellipses) into boxes which are aligned along the initiated 
bar. A particle in a box orbit is unable to precess all the 
way round and will fall each time back to the bar. This ef- 
fect will cause the bar to increase in both size and strength. 
To measure the radial stability of the systems we fitted the 
shape of an ellipsoi dal mass distribution by means of an 
iterative procedure dDubinski fc Carlberg 199ll ; iKatz 199ll ; 
iMeza fc Zamorano 1997 ; Meza 2002! ) at every half-mass dy- 



namical time. This detects any bar feature that is located 
within a given radius. The initial condition of this method 
is 
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0.95 




0.85 



Figure 4. Axis ratios as a function of time of the anisotropic 
systems with a decreasing anisotropy. 



y 2 z 2 ^ 1 ' 2 
p = p(a) with a = ( x 2 + — + 



(21) 



and with Mij = ~^ L i the principal components of the 
inertia tensor M zz ^ M yy ^ M xx and the axis ratios q and 
s equal to 1, assuming a spherical mass distribution within a 
certain sphere with a given radius for which we chose r — 5. 
For all considered models this radius encloses approximately 
70% of the total mass. To achieve these conditions a transi- 
tion to the center of mass has to be made followed by swap- 
ping the coordinate axes into the correct order. In the next 
step the eigenvalues and eigenvectors of the inertia tensor 
lij are calculated, transforming it into a diagonal matrix. 
At this point the new axis ratios can be calculated 



My_ 

M x 



1/2 



= — and 

a 



M x 



1/2 



(22) 



which in turn are used as the conditions for the next iteration 
step. The iteration was stopped as soon as both axis ratios 
converged to a value within a pre-established tolerance of 



10" 



Thus at each half-mass dynamical time the values of 



these axis ratios serve as measures of the strength of the bar 
instability, if present. 



3 HERNQUIST MODELS WITHOUT A 
BLACK HOLE 

In this section we investigate the stability of two different 
families of anisotropic Hernquist models without a central 
supermassive black hole. For the analytical construction of 
these models we refer to Baes & Dejonghe (2002), however 
we will recapitulate the characteristics of each family. 



3.1 Family I: Decreasing anisotropy 

We find Hernquist models with a decreasing anisotropy by 
assuming an augmented density of the form 
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Figure 5. The stability of the Hernquist models without a black 
hole and with increasing anisotropy (see Section 13.21 1, expressed 
as the minima of the axis ratios c/a during the simulations. The 
shaded area indicates the region of physical systems, i.e. with a 
non-negative distribution function. 



p(lp,r) 



4-2/3„ 



(1 + r) 



2tt (i - ipy-^o 



p 2/3 



(23) 



with (in = Po 



and n a natural number. After some 



algebra we find the distribution function 

2 /3o 



F(E,L) = 



(2tt) 5 /: 



T(5 - 2p n )L- 2f3a E 5/2 ~ 2 ^ +f3 ° 



E 



i 



oW r ( 2 -A,)r(2^-2/?„+/3 ) 

x 2 F 1 (5-2p n ,l-2p ;^-t-2p n +p ;E) ,(24) 



with 2-Fi a hypergeometric function (see Appendix A), and 
the anisotropy 



P(r) = 1 - 



aUr) _ Po + Pur 



(r) 



1 + r 



(25) 



which decreases as a function of radius. Since for Po ^ 
we only find tangentially dominated systems which are free 
of radial instabilities, we limit ourselves to the investigation 
of the case Po = 0.5. For this value, n = corresponds to 
a system with constant anisotropy. We plotted the axis ra- 
tios c/a for a number of different models with different n in 
Fig. [4] Here and in the remainder of the paper, we define 
those models that keep the axis ratio c/a > 0.95 over 50 
dynamical times as being stable. The only model that does 
not satisfy this criterion is that with n = 1, which is every- 
where radially anisotropic. For n ^ 2, the models become 
tangentially anisotropic at larger radii and as a consequence 
are much more stable. This is evident from Fig. [4] It is clear 
that the minimum of c/a is reached rapidly, whereafter the 
systems are in an equilibrium state, but are slightly non- 
spherical. 
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3.2 Family II: Increasing anisotropy 

These models are a gene ralization of the Osipkov-Merritt 
models i|Cuddeford 1991) with an augmented density and 
DF of the general form 



p(i[>,r) = r" 2,3o /W>) (1 + Ar 2 ) 1+00 with A = -^, (26) 



F{E,L) = F (Q)L-^° with < Q = £ - — < 1, (27) 

2r a 

and E denotes the energy, L the angular momentum and 
r a the anisotropy radius. The explicit form of f(ip) for 
the Hernquist potential-density pair can be found in Baes 
& Dejonghe (2002). As mentioned by them, the DFs can 
be written analytically for the half-integer values /3o = 
0.5,0,-0.5,-1, so we will limit ourselves to these cases. 
For every value of flo, we also computed numerically the 
maximum anisotropy value A max (/9o), outside which the DFs 
become negative for some values of Q and L. The area of 
physical systems is indicated in Fig. [S] Our models have the 
following functional form: 



• ft, = 0.5 

F(E,L) = 

• A) = o 

F(E,L) = 



Q 3Q 2 + A(3Q 2 - 5Q + 2) 
3 arcsin yfQ 



(28) 



8\/27r 3 



with 

HQ) 



• A> = 

F(£,L) 



+ v/Q(l-2Q) 



(i-O) 2 



+ 8A 



. A = -0.5 
F(S,L) 



Lf(Q) 



256\/27r 3 (l - Q) 5 

" HQ) 



■. arctan 



( ^ ) 


, HQ)] 


Wi-QJ 


VQ \ 



(29) 



(30) 



4tt 3 (i - Q)V<3 2 + A ( 1 -0) 2 ' 

6(1 + A) 2 Q 6 - 2(16A 2 + 26A + 10)Q 5 
+(70A 2 + 87A + 20)Q 4 - 2A(40A + 33)Q 3 
+A(50A + 19)Q 2 -16A 2 Q + 2A 2 . (31) 



(32) 
(33) 



with 

fi(Q) = 15 [(16A + 120)Q 2 - (72 + 32A)Q + 15 + 16A] , (34) 

HQ) = 384(1 + A) 2 Q 6 

-(1984A 2 + 3712A + 1728)Q 5 

\ 2 



(1-Q) 5 / 2 



+(4160X Z + 7008A + 2784)Q 4 
-(4480A 2 + 6192A + 1200)Q 3 
+(2560A 2 + 2368A + 930)Q 2 
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Figure 7. Contour plot of the axis ratio c/a for the Osipkov- 
Merritt models as a function of anisotropy radius and mass of the 
SBH. The shaded area indicates the region of physical systems, i.e. 
with a non-negative distribution function. The solid lines indicate 
the minimal values during the simulation, the dashed lines show 
the axis ratios at the end of the simulation (at t = 50T h ). 



models. In this respect a dark matter halo and a central 
supermassive black hole are important and can change the 
galaxy's properties dramatically. However, up to now there 
are few known analytical systems that include e.g. a super- 
massive black hole. The only models known so far are pre- 
sented in Ciotti (1996), Baes & Dejonghe (2004) and Baes 
et al. (2005) which are all based on the 7-models with spe- 
cial attention to the Hernquist model and in Stiavelli (1998) 
where the distribution function of a stellar system around 
an SBH is derived from statistical mechanic considerations. 

In this section we investigate the radial stability of both 
isotropic and anisotropic Hernquist models containing a su- 
permassive black hole, as these represent the closest analyt- 
ical approach to the observations. For the following sections 
we will use the representation of Ciotti (1996); again, we are 
not going into great detail in the derivation of the analytical 
distribution function. 

In essence the DFs are obtained from an analytical 
Osipkov-Merritt inversion of the systems governed by eq. ([2} 
and ([3]). As a consequence, these models can be viewed as 
a extension of eq. (|29|l . Subsequently, we will refer to these 
combined systems as Osipkov-Merritt models. The DFs can 
be written as 



F(Q) = Fi(Q) + 



F a (Q) 



(37) 



-(704A 2 + 240A + 225)Q + 64A 2 



(35) 



For all models the anisotropy is given by the simple formula 



/3(r) = 



■ 2 + fori 
r 2 + r\ 



(36) 



showing an increase in anisotropy as a function of radius. 
The results of the A-body investigation for all /3q and A are 
summarized in figure [S] where we plotted the minimal axis 
ratios c/a for the DFs in this parameter space. To derive 
this plot, we simulated systems with j3o = 0.5,0,-0.5,-1 
and A = 1,2,4,6,10,16,24 where physically possible. The 
case where /3o — corresponds to the traditional anisotropic 
Osipkov-Merritt Hernquist model that has been previously 
investigated in a similar way by Meza & Zamorano (1997). 
These authors state the system with r a w 1.1 (or A ~ 0.82) 
as stable. To compare our study with theirs, we simulated 
this model in addition to the other systems. For this model 
we find an axis ratio c/a « 0.95 after 50 dynamical times, 
and 2K r /K t ~ 2.2 during the entire run. These values are in 
agreement with their results, therefore we will define c/a = 
0.95 as our stability criterion. 

As is to be expected, the anisotropy radius A strongly 
affects the formation of radial-orbit instabilities, so that only 
models with a low value of A remain stable. Furthermore we 
note that all models remain in their new equilibrium state 
after t ~ 107^, as in case of the DFs of Family I. As an 
example, the c/a ratio evolution for one of the systems is 
given in Fig. [3] 



4 HERNQUIST MODELS WITH A 
SUPERMASSIVE BLACK HOLE 

The now established presence of diverse components in a 
great variety of galaxies calls for more advanced dynamical 



where Q has the same definition as in eq. 
ural parameter q is defined through 



= Q 1 



/' 



271, A more nat- 



(38) 



4.1 Family I: Isotropic 

We find an isotropic system by letting r a diverge to 00. Then 
eq. (|37[l simplifies to 



F(E) = E(E) 



2V8V 3 



dE 
H 



d 
Jl 



Ft (I) 



(39) 



with the argument / defined as I 2 = 1 — q. For F^Q,) we refer 
to Ciotti (1996) as this involves combinations of elliptic and 
Jacobian functions. These models only differ from those of 
Baes et al. (2005) in the definition of the parameter fi. 

Although the systems are isotropic, their DFs have a 
local maximum when /i > (as shown in Fig. [6jl . Hence, the 
sufficient criteria of Antonov (1962) and Doremus and Feix 
(1973) for isotropic systems cannot be applied. However, in 
our subsequent analysis of the systems with and without 
an SBH in Section 5, it will be shown that all models with 
r a > 1 are stable. In other words, it becomes evident that the 
addition of a central SBH, although it changes the dynamics 
dramatically, does not influence the density distribution of 
an isotropic system. 

4.2 Family II: Anisotropic 

In a similar way as the isotropic case the distribution func- 
tion can be written as 



F(Q) 



Fi(Q) - 
1 

2^871-3 



Fa(Q) 



dQ 
dl 



d_ 
~dl 



Ft(l) + ^fi 



(40) 
(41) 
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H=0,r a =0.25 - n= 0.01. r a = 0.25 - n= 0.03, r a = 0.25 - n= 0, r a = 0.5 - ]l= 0.01, r a = 0.5 - H= 0.03. r a = 0.5 -H=0,r a =1.0 - n= 0,01. r a = 1.0 - n= 0.03, r a = 1.0 
(1= 0.05, r a = 0.25 - \l= 0.07. r a = 0.25 - \i= 0.1. r a = 25 - \l= 0.05, r a = 0.5 - (i= 0.07, r a = 0.5 - (1= 0.1 , r a = 0.5 - |l= 0.05, r a = 1.0 _ (1= 0.07, r a = 1.0 - |1= 0.1, i a = 1.0 
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Figure 8. The evolution of systems with an anisotropy radius of respectively r a = 0.25,0.50, 1.0, but with different mass of the SBH 
(fi = 0,0.01,0.03,0.05,0.07,0.1). For each column we plot the evolution of the axes ratios (dashed lines) and Fourier coefficients (full 
lines) in the upper three rows. The bottom row contains the evolution of the 2K±/Ku ratio as a function of dynamical time (T^). 



where again F^(l) is defined in Ciotti (1996). In Fig. |S]we 
display systems with several values of fi and r a . Notice that 
for small values of r a the DFs have a local minimum. As 
a consequence, for every fi there exists a smallest possible 
r a , where this minimum becomes zero; smaller values of this 
boundary r a result in negative DFs, thus creating unphys- 
ical systems. For fi = 0, the minimal anisotropy radius is 
r a ~ 0.202; for /i = 0.1 the boundary becomes r a ~ 0.240. 
From the viewpoint of a stability analysis these systems are 
the most interesting. In the following section we will discuss 



their evolution in detail, comparing them with the models 
without an SBH (eq. p5]l). 



5 STABILITY ANALYSIS OF THE 
OSIPKOV-MERRITT MODELS 

To investigate any trend about the radial stability of 
these systems, we investigate the 2-parameter space (r a ,/x). 
In total, we performed 25 simulations, with r a = 
0.25, 0.50, 0.70, 0.85, 1.00 and fj, = 0.01, 0.03, 0.05, 0.07, 0.10. 
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Figure 9. The spatial density p and the velocity dispersions <r± and try, in the three principal planes, of an Osipkov-Merritt Hernquist 
system without an SBH, and with r a = 0.25. Dynamical times t = 0, t = 257)j and t = 507)j are displayed. 



This way we derived a grid of values of the c/a axis ra- 
tios, shown in Fig. [7] As stated before, the (p = 0)-axis 
corresponds to the systems in eq. ((29}. The solid lines indi- 
cate the minimal values reached during the simulation (i.e. 
when the instability is strongest). The rate at which these 
minima are reached strongly depends on r a , ranging from 
a few half-mass dynamical times for highly radial models 
to t ~ 50Th for systems with r a = 1.00. After the point of 
time upon which a system obtains its minimal c/a the influ- 
ence of the SBH causes a diminution of the bar instability, 
resulting in the c/a axis ratios at t = 507), shown by the 
dashed lines. Thus, in each system the particles are affected 



by two counteracting forces: the (relatively fast) bar forma- 
tion and the (more gradually) scattering near the center due 
to the spherically symmetric gravitational potential of the 
black hole. The contour line c/a — 0.95 is highlighted as our 
stability criterion. 

A full dynamical analysis would require a detailed study 
of the orbital distribution of the stellar mass. However, we 
can gain important insights into the dynamics of the models 
by visualizing the evolution of various quantities, in Figs.JSJ- 

H 

First, in order to retain a notion of 'radial' and 'tan- 
gential' motion in an evolved system (resembling a triaxial 
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Figure 10. The spatial density p and the velocity dispersions a± and cry, in the three principal planes, of an Osipkov-Merritt Hernquist 
system with an SBH of p = 0.05, and with r a = 0.25. Dynamical times t = 0, t = 257), and t = 507), are displayed. 



model) at a certain time t , we use the method described in 
Section \2. 41 to approximate the mass distribution inside the 
radius n of each particle by an ellipsoid. Then, the velocity 
of a particle can be written into two components perpendicu- 
lar resp. parallel to its surface Vi — Vi t ± • Subsequently, 
the perpendicular and parallel velocity dispersion of the m 
nearest neighbors around a position r are 



_(r) = 



— E(v-«i)' 



(42) 



a " (r) = ^3T)E(^ii--ii)' 

In a similar manner we define 



2 

i,_L> 



i=l 

1 r 



(43) 



(44) 



(45) 



so that 2K±/K^ can serve as a non-spherical extension of 
2K r /K t . 
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Figure 11. The spatial density p and the velocity dispersions cr± and cry, in the three principal planes, of an Osipkov-Merritt Hernquist 
system with an SBH of p = 0.1, and with r a = 0.25. Dynamical times t = 0, t = 257^ and t = 507), are displayed. 



In Fig.[8]we compared our c/a criterion to other known 
methods to quantify the stability of the system. For all 
our systems with r a = 0.25, 0.50 and 1.0, we plot the 
evolution of the three axis ratios against the Fourier co- 
efficient (S ellwood fc Merritt 1994h in the three principal 
planes (XY, XZ and YZ), and the evolution of 2K±/K\\ (t). 
The Fourier coefficient in a plane is defined as 

3 

with 0j the azimuth of particle j in this plane. Together 
with the 2K±/Ku ratio it is widely used as a measurement 



of the instability of a system. From this figure it is clear 
that all 3 methods indicate a similar result. For the models 
r a = 1.0 (right column) the axis ratios and the 2K ±/K ii (t) 
ratio remain nearly constant during the entire run, while 
the Fourier coefficient either remains constant or increases 
upto a very low value and then slowly declines, indicating no 
further evolution in the stability has to be expected. Only 
the systems with fj, = and fi = 0.01 have moderately 
increasing Fourier coefficients, indicating a marginal, slowly 
evolving instability; models with higher fj, can be considered 
stable. 

Lower values of the anisotropy radius r a leads to un- 
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Figure 12. The spatial density p and the velocity dispersions a± and cry, in the three principal planes, of an Osipkov-Merritt Hernquist 
system with an SBH of p = 0.05, and with r a = 0.50. Dynamical times t = 0, t = 257), and t = 507), are displayed. 



stable systems; in columns 1 and 2 we display all systems 
with respectively r a = 0.25 and r a = 0.5. Clearly, the rate 
at which the bar is created depends strongly on r a . Dur- 
ing a simulation of an unstable system the 2K± / K\\ (t) ratio 
rapidly reaches a maximum and then declines to a value of 
ss 2.0 after which it remains constant for the remainder of 
the run. This maximum is an artifact of the softening applied 
in the N-body simulation and is caused by the fact that the 
initial conditions are not exactly in dynamical equilibrium 
(see Sect. 2.3). Lowering the time-step and softening length 
caused a diminution of the maxima. The time at which the 
2K±/K\\ ratio reaches this constant value corresponds with 



the time the axis ratios reach their lowest value and the 
Fourier coefficients reach their highest value. The influence 
of the mass of an SBH is more clearly visible in the axis 
ratios and Fourier coefficients. A remarkable result is the 
difference between the systems r a = 0.25 and r a = 0.5 in 
the evolution of the axis ratios: the former models first be- 
come triaxial, after which b/a increases and c/b decreases. 
The latter models are prolate axisymmetric during the entire 
run. 

In summary, an SBH mass of a few percent can 
prevent or reduce the bar instabilities in anisotropic 
systems. This result agrees well with similar studies in 
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disk galaxies dNorman et al. 19961; IShen fc Sellwood 2004 



lAthanassoula et al. 20051 ; iHozumi fc Hernquist 20051 ). The 
effect is strongest for models with strong radial anisotropies, 
where the decrease of the bar strength is proportional to the 
SBH mass. In other words, while more radially anisotropic 
systems develop stronger bars than more isotropic models, 
the bars of the former are more easily affected by a super- 
massive black hole (see Fig. [8}. This is to be expected, since 
radial systems host more eccentric orbits, therefore more 
particles from the outer regions pass near the center where 
their orbits can be altered by the Kepler force of the SBH. 

An alternative approach to present the evolution of the 
models is given in Figs. 1 91121 where we show the evolu- 
tion of 4 systems by means of the density p(r) and ve- 
locity dispersions <Tj_(r) and cry (r), in at dynamical times 
t = 0, t = 25Th and t = 502^. In each principal plane the 
moments are calculated on a grid of 2500 points, with 50 
nearest neighbors around every grid position. 

Fig. [9] displays an Osipkov-Merritt system without an 
SBH and anisotropy radius r a = 0.25. As was shown in Fig. [8] 
this model has a strong bar formation, resulting into a new 
equilibrium state after t — 5Th which it retains during the 
rest of the run (as can be seen at t = 25Th and t = 507^). 
This bar alters the density distribution into a roughly triax- 
ial symmetry, even peanut-shaped in the XZ-plane where 
the radial instability is the most prominent. The tangential 
dispersion an increases significantly. This occurs especially 
at the edges, where in contrast the radial dispersion van- 
ishes. This can be explained by the mechanism of the bar 
formation: particles that pass through the bar are pulled to- 
wards it, and eventually align their orbit with the bar. Only 
the orbits along the principal axes remain largely unaffected 
by the bar due to the symmetric forces on these particles, 
hence their motion remains radial. 

In Fig. [TO] a model with r a = 0.25 and fi — 0.05 is 
shown. Again a bar is formed, but less pronounced than in 
the absence of an SBH. Clearly, during the run the bar is re- 
duced by the SBH, causing a gradual increase in the c/a axis 
ratio (XZ-plane). More striking however is the evolution in 
the XY-plane, where the ellipticity has disappeared. Thus, 
the model has become an oblate axisymmetric system. This 
is also reflected in the dispersions: en again follows the bar 
structure, but the cross-form a± vanishes as particles pass 
near the SBH. Since most particles reside in the XY-plane, 
on eccentric orbits (since r a is small), the scattering in this 
plane is strongest. After t = 507^, we expect a further small 
increase in the c/a axis ratio, but as the velocity dispersion 
becomes more isotropic fewer particles from the outer re- 
gions will pass near the center (i.e. be affected by the SBH), 
hence the model will not change much further. 

This can also be seen by comparing the system with 
/i = 0.05 to a model with /i = 0.1 (Fig. Ill[) . This model has 
essentially the same properties as the former. The larger 
SBH mass has above all influence on its efficiency, resulting 
in a faster bar reduction. 

Finally, we consider the effect of the anisotropy radius 
by analyzing a system with [i — 0.05 and r a = 0.5 (Fig. ll2[) . 
Compared to Fig. 1101 the initial bar is less strong, as ex- 
pected. However, its structure and evolution is different from 
the system with r a = 0.25. First, a± remains spherically 
distributed during the run, thus less scattering occurs. This 



implies less reduction of the bar instability. Moreover, the 
density does not become symmetric around the Z-axis. In 
contrast, the X-axis is now the symmetry axis during the 
entire run, resulting in a prolate axisymmetric system. It 
thus seems that models with an SBH become oblate or pro- 
late, depending on their velocity anisotropy. It would indeed 
be very interesting to compare the orbital structure of both 
these systems in full detail. 

As a final remark we note that inside a radius tk = 
y / /I/(l — \fP) the force of the SBH is stronger than the stel- 
lar component, so that all models remain spherical inside 
this radius. In conclusion, systems with an SBH become ax- 
isymmetric systems with a spherically symmetric core. 



6 CONCLUSIONS AND SUMMARY 

Most mass estimates of SBHs result from dynamical mod- 
els of either stellar or gas kinematics. The inclusion of 
strong radial anisotrop y is considered in these models 
(jBinnev fc Mamon 19821 ) . yet they have never been tested 
for radial stability. Our goal was to test the stability of 
systems with a central SBH and to look for any trend 
as a function of the mass of the SBH. We used the 
same method that was previously introduced by Meza & 
Zamorano (1997) and extended it to systems with a central 
SBH. We first tested th e procedure on Hernquist systems 
jBaes fc Deionghe 2 002) without an SBH and with different 
anisotropic behavior. Our method appeared to be efficient 
in discriminating the stable from the unstable systems. 

Instead of focusing on complicated numerically derived 
dynamical models, we opted for analytical distribution 
functions that take the effect of a central SBH into account, 
in order to be able to look for any trend. Since the isotropic 
Hernquist models with an SBH do not have distribution 
functions that are monotonically increasing functions o f 
the binding energy l|Ciotti 19961 ; iBaes fc Deionghe 20041 ) 
and hence the sufficient criteria of Antonov (1962) and 
Doremus and Feix (1973) for isotropic systems cannot 
be applied, we first investigated the radial stability of 
these systems. No effect was found by letting the mass 
of the SBH vary, giving only stable systems. However, 
in the case of the anisotropic systems with an SBH we 
did find a dependence of the stability of the system on 
the mass of the SBH. The more massive the SBH, the 
more stable a system becomes, but especially the more the 
instability is reduced. A trend which is most obvious in 
very anisotropic systems (thus with very small anisotropy 
radius r a )- An SBH with a mass of a few percent of the 
entire galaxy mass, is able to weaken the strength of the 
bar, which is in correspondence with similar studies in 
disk galaxies dNorman et al. 19961; IShen fc Sellwood 20041 ; 



lAthanassoula et al. 20051 Hozumi fc Hernquist 200 



Judging from Fig. [7] the stability boundary of c/a > 0.95 
over 50 dynamical times, shifts from r a ~ 1.1 for /i = to 
r a w 1.0 for fi = 0.1. This corresponds to 2K r /K t = 2.2 for 
H — and to 2K r /K t = 2.0 for /i = 0.1. These values are 
in very good agreement with previous authors. 

Remarkably, systems with an SBH but with different 
anisotropy radii r a evolve differently: highly radial systems 
first become triaxial whereafter the SBH makes them more 
oblate if fi is large, while less radial models tend to form 
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first into axisymmetric prolate structures, that then become 
less elongated due to the influence of the SBH. 

It is also interesting to note that the central density dis- 
tribution of systems with an SBH remains spherically sym- 
metric during the entire simulation out to a radius of half 
the effective radius. This is not the case for systems without 
an SBH, which become axisymmetric or triaxial, depending 
on r a . Interestingly, this includes the region that is con- 
sidered for the Mbh — a relation, which predicts such an 
evolutionary link between the central SBH and the spheroid 
where it resides. Similarly, the central anisotropy parame- 
ter decreases as a function of time at a rate proportional to 
the mass of the SBH, due to more tangential orbits at the 
center. 

Apart from a central SBH, one can also 
investigate the effec t of central density cusps 
l|Sellwood fe Evans 200ll; iHollev-Bockelmann et al. 200lT ) 
or isotropic cores |Trenti fc Bertin 20061 ? on the radial 



stability. Previous research shows that both central density 
cusps and isotropic cores act as dynamical stabilizers, hence 
the same effect as the addition of an SBH. In the future 
we plan to look at the combination of both investigations, 
namely the combination of an SBH and central density 
cusp. 

Ultimately, one can of course verify the radial stability 
of the state-of-the-art models that are being used to estimate 
the mass of the SBHs with e.g. our methodology, however this 
investigation was not the goal of this paper. 
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APPENDIX A: 
FUNCTION 



HYPERGEOMETRIC 



A relevant definition of Generalized hypergeometric func- 
tions p F q (ai, . . . a p ;bi, . . . ,b q ; x) can be found in Gradshteyn 
& Ryzhik Sec. 9.14, page 1071 in the 5th edition. For specific 
arguments hypergeometric functions can be reduced to more 
simpler analytical functions, this can be done with mathe- 
matical software packages as e.g. Maple or Mathematica. 
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For general coefficients however, the evaluation needs to be 
done by means of hypergeometric series: 

+00 

p F q (a 1 ,...,a p ;b 1 ,...,b q ;x) = 1 + ^ ( A1 ) 

3=0 

n xU P i=l( a i + k ) 
(!=0 (i+*)n?=i(6i+*)' 

This routine is very suitable to numerically evaluate any 
generalized hypergeometric function within a certain degree 
of accuracy. 



